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Abstract 

We report a novel approach for inversion of large random matrices in massive Multiple-Input Multiple Output (MIMO) 
systems. It is based on the concept of inverse vectors in which an inverse vector is defined for each column of the principal 
matrix. Such an inverse vector has to satisfy two constraints. Firstly, it has to be in the null-space of all the remaining 
columns. We call it the null-space problem. Secondly, it has to form a projection of value equal to one in the direction of 
selected column. We term it as the normalization problem. The process essentially decomposes the inversion problem and 
distributes it over columns. Each column can be thought of as a node in the network or a particle in a swarm seeking its 
own solution, the inverse vector, which lightens the computational load on it. Another benefit of this approach is its 
applicability to all three cases pertaining to a linear system: the fully-determined, the over-determined, and the under- 
determined case. It eliminates the need of forming the generalized inverse for the last two cases by providing a new way to 
solve the least squares problem and the Moore and Penrose's pseudoinverse problem. The approach makes no assumption 
regarding the size, structure or sparsity of the matrix. This makes it fully applicable to much in vogue large random matrices 
arising in massive MIMO systems. Also, the null-space problem opens the door for a plethora of methods available in 
literature for null-space computation to enter the realm of matrix inversion. There is even a flexibility of finding an exact or 
approximate inverse depending on the null-space method employed. We employ the Householder's null-space method for 
exact solution and present a complete exposition of the new approach. A detailed comparison with well-established matrix 
inversion methods in literature is also given. 



Citation: Anjum MAR, Ahmed MM (2014) A New Approach for Inversion of Large Random Matrices in Massive MIMO Systems. PLoS ONE 9(4): e94958. 
doi:10.1371/journal.pone.0094958 

Editor: Shyamal D Peddada, National Institute of Environmental and Health Sciences, United States of America 
Received November 26, 2013; Accepted March 21, 2014; Published April 14, 2014 

Copyright: © 2014 Anjum, Ahmed. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Funding: The authors have no support or funding to report. 

Competing Interests: The authors have declared that no competing interests exist. 
* E-mail: ali.raza.anjum@gmail.com 



Introduction 

Multiple-Input Multiple-Output (MIMO) systems form a well 
established area of wireless communications [1]. A significant 
increase in interest in MIMO systems has been witnessed in the 
last few years due to the advent of massive MIMO systems [2,3]. 
In these systems, more degrees of freedom in terms of data rate 
and link reliability are available due to the increase in the number 
of transmit and receive antennas. These advantages become even 
more impressive in multi-user scenario when there is a possibility 
to transmit to several users simultaneously [2] . However, there is 
still a challenge of low complexity detection techniques for the 
practical realization of such systems [4] . Many high performance 
detection methods for massive MIMO systems require an 
unconstrained solution to a linear estimation problem [5]. Linear 
estimation requires the inversion of channel matrix which, in such 
systems, can be problematic because of its potentially large size. 
For example, matrices of size 40*40 and above have recendy been 
reported in literature as massive [2]. 

Therefore, there is a need to find solutions that do not require 
outright inversion. Various methods in this regard have been 
proposed in literature: Cayley-Hamilton method, Neumann series 
method, QR method, random matrix methods, LSMR, LSQR, 
Kaczmarz method and the ones based on polynomial and 
truncated polynomial filters [6-19]. These methods still require 
a lot of computational effort and some of them have even proven 



to be more complex than the outright inversion [2,6]. While it is 
possible that some may perform better than other, it is not always 
possible to have a fair comparison among them since these 
methods follow independent approaches and are dependent on 
different sets of parameters. 

Keeping that in view, our objective in this paper is to develop an 
alternate method to find the inverse of the channel matrix. While 
addressing this problem, we dispense with three fundamental 
assumptions that are generally implied in the traditional methods. 
We endeavor to do so in order to provide more room for thought. 
First, no assumption has been made about the structure of the 
channel matrix. Secondly, the matrix can be purely random 
instead of being deterministic. Finally, we do not assume that it is 
sparse [20] . Keeping this in view, we report a novel method that 
not only finds the inverse of a channel matrix but also brings a new 
perspective to the inversion process itself. The proposed method is 
a comprehensive one and is fully applicable to all three matrix 
inversion cases pertaining to a linear system. We also provide its 
detailed comparison with well-known matrix inversion methods 
available in literature for the sake of analysis and, hence, 
understanding. 

Rest of the article is organized as follows. To begin with, the 
subsequent section lays out basic system model and the essential 
nomenclature. The third section presents the basic idea behind the 
proposed method. The inversion problem is solved according to 
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the proposed method using Householder's null-space method in 
the fourth section. In fifth section, solution to three fundamental 
cases of linear systems is presented. Proposed method is then 
compared with the QR decomposition (QRD) method, Least 
Squares (LS) method, and Moore and Penrose's pseudoinverse 
method in sections six, seven, and eight respectively. A complete 
algorithm for step by step computation of the inverse matrix 
according to the proposed method is outlined in section nine. 
Simulation results demonstrating the speed and accuracy of the 
proposed method in comparison to the state of the art methods 
available in literature are presented in section ten. Finally, the 
article concludes with a brief discussion in section eleven. 

System Model 

A MIMO system is represented by a system of linear equations 
[21]. 

Hx = y (1) 

H is the channel matrix of dimension m x n, x is the transmitted 
vector of dimension n x 1 , and y is the received vector of 
dimension mxl. Entries of H matrix are complex-valued 
independent and identically distributed Gaussian random vari- 
ables of zero mean and unity variance [6]. m and n are the 
number of transmit and receive antennas respectively. If m = n, 
the number of equations is equal to the number of unknowns and 
the system is fully-determined. If m>n, the number of equations is 
greater than the number of unknowns and the system is over- 
determined. If m<n, the number of equations is less than the 
number of unknowns and the system is under-determined. 

Basic Idea 

We begin with the assumption that H has a full-rank. Re- 
writing Eq. (1) as: 

y = hjXj+ Y™i = \hiXi (2) 

First term on the right hand side of Eq. (2) denotes y'-th column 
of matrix H.Xj is the associated unknown. The second term 
represents the remaining columns of H with their respective 
unknowns. This term can be separately written in matrix form by 
defining a new column reduced matrix Fj and a symbol reduced 
vector x. 

y = hjXj+FjX (3) 

Fj has dimensions of m x (n— 1). Taking the dot product of Eq. 
(3) with an arbitrary column vector Wj of dimension mxl, 

wfy = wfhj Xj + wfFjX (4) 



wfFj = 0 (6) 

Eqs. (5) and (6) define the j-th row of the inverse matrix. It 
should have a dot product of one with the j-th column and zero 
with all the remaining columns. The same is also true for all the 
other rows of the inverse matrix. Once we are able to solve Eqs. (5) 
and (6) for Wj, all rows of inverse matrix can be computed one by 
one in the same fashion and a complete inverse matrix can be 
built. Eq. (5) essentially restricts the length of the null-space 
solution produced by Eq. (6) and serves as a constraint for Eq. (6). 
Therefore, we proceed by solving the Eq. (6) first. It can be solved 
by a variety of methods available in literature of for computing the 
null-space of a matrix. We propose to solve it by Householder's 
null-space method because it is stable and provides an exact 
solution. 

Solution Using Householder's Null-Space Method 

Householder's method has been traditionally viewed as a means 
to achieve the QR decomposition (QRD). It performs a series of 
orthogonal transformations on an arbitrary matrix to convert it to 
an upper triangular matrix. These transformations can be 
performed either by Given's rotation matrices or Householder's 
reflection matrices. A reflection matrix can perform the job of 
many rotation matrices at once. Therefore, reflection matrices are 
more desirable in our case. We proceed by performing a series of 
orthogonal transformations on the matrix Fj using the reflection 
matrices to convert it to an upper triangular matrix R. n—\ 
reflection matrices will be required for that purpose because Fj has 
n — 1 columns. 



H„-i . . . H3H2H1F/ = QFj = R (7) 



Hj's represent the Householder's reflection matrices. They are 
square, symmetric and orthogonal with each one having the 
dimensions of m x m. Q has also the same dimensions. But the 
dimensions of R are mx(n— 1). This is because the first n— 1 
rows of Q generate an (n — 1) x (n — 1) upper triangular portion in 
R. Remaining m — (n— 1) rows in R are zero. These rows are 
produced by last m — (n — 1) rows of Q. 



R=[n r 2 ... r„_i 0 ... 0 0 } H 
Q = [ qi q 2 ••• q„-\ q„ ■■■ q m -\ q m ]" (8) 
R=[n r 2 ... r„_i 0 ... 0 Of (9) 



To evaluate the j-th unknown Xj from Eq. (4), 



Wfhj=\ 



(5) 



q\ mark the rows that produce zeros in R matrix. Substituting 
Eqs. (8) and (9) in Eq. (7), 
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Concentrating on ^'s only in Eq. (10), 
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q's form the basis for the left null-space of Fj and, hence, solve Eq 
(6). Only equation that remains to be solved now is Eq. (5). In 
order to do so, we would need to consider the three fundamental 
cases of linear systems. 

Three Fundamental Cases 

In this section, we take up the three fundamental cases of linear 
systems and solve them one by one. 

1. The Fully-Determined Case 

In the fully-determined case, there is only one vector in the left 
null-space of Fj, i.e., m — (n— 1) = 1 given that m = n. The solution 
in Eq. (1 1) is unique. 



(12) 



The only remaining task is to rescale qj according to Eq. (5). 

(13) 
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The identity matrix will have the dimensions of n x n. 

2. The Over-Determined Case 

In the over-determined case, the matrix H will already have a 
right null-space due to a greater number of rows than the columns. 
The removal of extra column to form Fj will populate the null- 
space even further and there will be m — (n— 1) vectors in the left 
null-space of Fj. Hence, the solution to the Eq. (6) will not be 
unique. 



kj = qfhj 



(14) 



kj is the rescaling factor for the j-th row of inverse matrix Fj. By 
iterating j for all columns of H matrix, complete inverse matrix W 
can be built. 
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All the vectors from q n to q m solve Eq. (6). An obvious question 
would be which one to choose. At this point, we adjourn this 
question to Section VIII where it would be amply addressed. For 
now, any one of them can be selected. The next step would be to 
rescale it according to Eq. (13). By repeating the same procedure 
for all columns of//, a complete inverse matrix can be built. 

3. The Under-Determined Case 

In the under-determined case, H matrix will have right null- 
space due to a greater number of columns than the rows. In order 
to be consistent with the earlier work, we can transpose it and 
move the null-space to left. The inverse matrix can then be 
computed in the manner described in subsection (b). Only a 
transposition of inverse matrix is required afterwards to revert to 
the original case. 

Comparison with QR Decomposition 

Traditional methods employing QR decomposition take a 
linear system of the form, 



into, 



Ax = b 



x = R- l Q H b 



(17) 



(18) 



by factorizing A in Eq. (17) into an orthogonal matrix Q and an 
upper- triangular matrix R [22]. There are two famous ways for 
achieving QR decomposition: Gram-Schmidt QR and House- 
holder's QR [22]. Traditionally, Gram-Schmidt method has been 
used for that purpose. Objective was to convert anmx/i matrix A 
into an m x n orthogonal matrix Q. Matrix R connects A to Q. 



(19) 



Householder's method takes a shift in perspective. R is targeted 
instead of Q. Being an upper triangular matrix, R is much more 
suitable for back substitution. Q serves as the connecting matrix. 



(20) 



The Q in Eq. (20) has extra m — n columns than the one 
obtained by Gram-Schmidt method. First m columns serve as the 
orthonormal basis for the column space of A. Remaining m — n 
columns form the basis for the left null-space of A because R can 
only be made triangular upto first n rows. [22]. 



A mn — [ Qmn Qm{ni — n) 



R„ 



{m — n)n 



(21) 



In the proposed method, Househoder's method is employed to 
find the null space basis Q m ( m -n) °f Fj. Also, there is a difference 
in the orientation of Q in Eq. (8) which leads to, 

A = Q H R (22) 



instead of A = QR. We take Q as such because its last m — n 
rows are responsible for the zero last m — n rows of R. 
Orthonormal column space of Q in Eq. (20) is now an 
orthonormal row space in Eq. (22). 

Comparison with Least Squares 

When A has fewer columns than the rows, the system becomes 
over-determined. There are more equations than the unknowns. 
The solution is then obtained by Least Squares (LS). 



A H Ax=A H b 



(23) 



QR decomposition is a popular method to solve LS problems 
and when applied to Eq. (23) will result in Eq. (18). 

We have not approached the over-determined case using LS. 
We have taken it column by column. A dearth of columns in an 
over-determined matrix signals an apriori presence of the left null- 
space. Removal of one more column will populate it further. This 
will bring us to select, out of this abundance of null-space vectors, 
the one which has the highest value of A: depicted in Eq. (14). It 
needs to be so because the selected vector will then be divided by k 
according to Eq. (13). Smaller values of k can inflate the 
magnitude of the null-space vector after division. In the worst 
case of k = 0, there will be a division by zero. Therefore, out of 
many null-space vectors obtained by Householder's null-space 
method, the best would be the one which has the highest value 
of k for the specific column vector. The term best can be pushed 
even further. In ideal case, a value of one for k would imply no 
division at all. But the caveat is to avoid smaller values of k in 
order to prevent instability. 

Comparison with Moore and Penrose's Pseudoinverse 

Now we turn to the last case, the case when the matrix A is 
singular. A solution in this case is sought by forming the Moore 
and Penrose's pseudoinverse. Not an exact inverse though, the 
idea is to find the shortest vector that solves the LS equation. The 
term shortest implies a vector that has no null-space component. 
One way to do is to use Singular Value Decomposition (SVD) 
[20]. SVD decomposes the matrix A into two square, orthonormal 
matrices U and V. These matrices contain the basis for the 
column and row space of A respectively. 



A=UJ2 v H 



(24) 



The third matrix is a diagonal matrix that contains the 
singular values. Some of the singular values in w iU De zero f° r 
this case due to the presence of dependent columns in A. Columns 
of U and V for the corresponding zero singular values will 
represent the left and right null-spaces of A. Pseudoinverse is 
formed by leaving the zero singular values unchanged while 
inverting the others. This gives the shortest solution in row space. 

At this point, a 2 x 2 matrix with all entries equal to one would 
serve best to explain our approach. Once the first column is 
removed, Eq. (6) will seek a solution in the left null-space of the 
remaining columns. Eq. (5) will dictate this solution to have a 
projection of value equal to one with the removed column. This is 
impossible because both the removed column and the remaining 
column are exactly the same. A vector cannot be orthogonal and 
parallel to itself at the same time. In our opinion, this is precisely 
the reason for non-invertibility of a matrix rather than the 
traditional zero-determinant view. The only unique vector is the 
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Table 1. Comparison of the proposed algorithm with the state of the art methods available in literature. 
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vector left in the column space of F, the only remaining column. 
Therefore, it is rescaled with itself to satisfy Eq. (13). This is exactly 
the solution obtained by applying the SVD to above problem. 

Algorithm 

We summarize the steps discussed in the form of a complete 
algorithm as follows. 

1. If m>n, continue to step 2, otherwise transpose the matrix H 
before moving to step 2. 

2. Locate the 7-th column whose inverse vector is to be 
determined. 

3. Remove that column from H matrix to form a new matrix Ft, 

4. Decompose Fj to form Q and R matrices according to Eq. (7). 

5. Discard the R matrix and the first n rows of Q matrix. 

6. Remaining m — n rows will constitute null-space of Fj. In case 
m = n, there will be only one row and, hence, one candidate for 
the inverse channel vector. In case m > n, all of the remaining 
m — n rows satisfy the criteria in Eq. (6). In this case, pick any 
one at random. 

7. Normalize the selected vector according to Eqs. (13) and (14) to 
form the inverse vector w?. 

8. Go back to step 2 and repeat the same procedure for the 
remaining columns of H matrix. 

9. Stack the inverse channel vectors on top of each other in the 
manner described in Eq. (15) to form the complete inverse 
matrix W . 

10. For the case when m<n, transposition of W formed in step 
9 is necessary to form the right handed inverse. 

A user-friendly code of this algorithm is provided in Code SI. 

Simulation Results 

We now present the simulation results of the algorithm. The 
components of the channel matrix are chosen to be independent 
and identically distributed (IID) circularly symmetric Gaussian 
random variables with a zero mean and unity variance, m is 
selected to be equal to n as this refers to multi-user case in M- 
MIMO systems because both the number of transmitting and 
receiving antennas become very large. For example, a typical 



matrix size can be 40 x 40 and above [2] . Also when m = n, an 
exact solution is possible and the residue and hence the error in the 
estimate can be zero. 

For simulation purpose, various matrix sizes have been selected 
and the results obtained in terms of the norm of residue, norm of 
the error in estimate and the computational time taken are 
displayed in Table 1 against the state of the art algorithms 
available in literature: LSMR, LSQR, and Kaczmarz method. 
Since the proposed method achieves zero error/ residue norms, its 
respective columns are not displayed in the table. Codes required 
for simulation of LSQR and LSMR algorithms have been 
downloaded from the website of Stanford University's System 
Optimization Laboratory and are used as is. 

Simulation results demonstrate that when LSMR, LSQR, and 
Kaczmarz method are applied to this problem, these methods 
suffer from very large residue norms. But that is not all. The norms 
for the error in estimate are even worse. Kaczmarz method is not 
only the slowest of all but also has largest error/ residue norms. 
Hence, the estimates become practically useless. On the other 
hand, the proposed method achieves perfectly zero error/residue 
norms in moderate time. Therefore, the proposed method can be 
a much better choice in this context due to its better accuracy and 
speed. 

Conclusion 

A novel approach for matrix inversion has been presented in 
this paper. The matrix was split into a set of column reduced 
matrices, each with its own inverse channel vector Wj. In order to 
qualify for an inverse channel vector, Wj had to satisfy two 
constraints; Wj I Fj = 0 and Wjhj = l. Householder's null-space 
method was used to solve the first constraint and scaling was used 
to solve the second one. Inverse matrix was then built from these 
inverse channel vectors. A detailed comparison with the state of 
the art methods available in literature was carried out all along to 
emphasize the novelty of the method. 

Supporting Information 

Code SI A user-friendly code of the proposed algo- 
rithm. 

(ZIP) 
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